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We improvise a novel approach to carry out first-principles simulations of graphene-based verti¬ 
cal field effect tunneling transistors that consist of a graphene|/i-BN|graphene multilayer structure. 
Within the density functional theory framework, we exploit the effective screening medium (ESM) 
method to properly treat boundary conditions for electrostatic potentials and investigate the effect 
of gate voltage. The distribution of free carriers and the band structure of both top and bottom 
graphene layers are calculated self-consistently. The dielectric properties of /i-BN thin films sand¬ 
wiched between graphene layers are computed layer-by-layer following the theory of microscopic 
permittivity. We find that the permittivities of BN layers are very close to that of crystalline /i-BN. 
The effect of interface with graphene on the dielectric properties of /i-BN is weak, according to an 
analysis on the interface charge redistribution. 


I. INTRODUCTION 

Graphene-based field-effect transistors have been 
synthesized,^ but the ratio of high and low resistances 
(switching ratio) is limited (less than a factor of 10) be¬ 
cause of the zero-gap band structure of pristine graphene. 
A possible solution is to introduce an energy gap into 
graphene by using, for instance, bilayer graphene, 
nanoribbons,^ or chemical derivatives.® Recently, Brit- 
nell et reported an alternative transistor architec¬ 

ture, a vertical field effect tunneling transistor (FETT), 
exhibiting a fairly high (^ 50) switching ratio. In FETT 
devices, a graphene|barrier|graphene trilayer is utilized 
as the current-carrying channel, as in Fig. 1(a); the elec¬ 
tric current flows into the channel through one graphene 
layer and out through the other. Inside the trilayer, elec¬ 
trons cross the thin barrier via tunneling. The chemical 
potentials of the two graphene layers and the tunneling 
conductance can be tuned by gate voltages. The bottom 
graphene layer (closer to gate electrode, Gre) can only 
partially screen the gate voltage because of the low den¬ 
sity of states (DOS) near the Dirac point, and so the top 
graphene layer (away from gate electrode, Grx) is also 
affected by the gate voltage. 

In a typical field-effect device, the current-carrying 
channel is separated by a thick dielectric from a metallic 
gate electrode. A gate voltage applied between the gate 
electrode and the current-carrying channel controls the 
free carrier density in the channel. The spatial distribu¬ 
tion of free carriers can be obtained by solving the electro¬ 
static Poisson equation. Macroscopic physical quantities 
of device components are used as parameters for Pois¬ 
son’s equation, including dielectric constants and elec¬ 
tron affinities of dielectrics and work functions of gate 
electrodes, etc.® However, the extrapolation of the macro¬ 
scopic theory of electrostatics to the nanometer scale 
is unjustified, because interfaces often exhibit dielectric 
screening properties different from the bulk.^® The ap¬ 
plicability of the macroscopic theory should be examined 
by studying the materials-specific dielectric screening 
properties of interfaces, one of motivations of this work. 
The density functional theory (DFT) method fully takes 


the atomic structure of a device into account, and mi¬ 
croscopic electrostatics quantities such as potential and 
charge density can be solved self-consistently, from which 
microscopic dielectric properties can be deduced. One 
can extract an effective dielectric constant from micro¬ 
scopic calculations to represent the dielectric screening 
of the interface, which is not necessarily a constant but 
a quantity that depends on the thickness of the interface 
and other factors. 

The boundary condition is crucial in solving for the 
electrostatic potential. In an FET device, for exam¬ 
ple, the electrostatic potential at the surface of a metal¬ 
lic gate electrode should be constant, and the electric 
field (derivative of the potential) in vacuum far away 
from the device should be zero. This mixed boundary 
condition is different from the conventional and widely- 
implemented periodic boundary condition. Otani et 
al}'^ proposed a Green’s-function-based effective screen¬ 
ing medium (ESM) approach to solve the electrostatic 
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FIG. 1. (Color Online) (a) Schematic of a graphene- 
based field effect tunneling transistor. The core graphene]/i- 
BN|graphene structure is separated from the doped Si (which 
serves as the gate electrode) by /i-BN crystal and SiOa, slabs 
(which serve as dielectrics). The graphene layers closer to 
and further away from the gate electrode are labelled as Gre 
and Grx, respectively, (b) The model used in our simulations. 
The dielectrics (h-BN and SiO^;) are replaced by vacuum. The 
core graphene I/i-BN I graphene structure is embedded betwen a 
semi-infinite vacuum {z > zi) and an ideally metallic medium 
{z < —zi, which replaces the doped Si gate). 
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potential under several different boundary conditions, 
which is promising for simulating the electric-field effect 
in planar devices. 

In this work, graphene-based FETT devices were sim¬ 
ulated using a DFT-I-ESM method. The computational 
approach, based on first principles, enables us to under¬ 
stand interface effects quantitatively, and thus enables 
computational design of functional systems. The rest 
of the paper is organized as follows. In the next Sec¬ 
tion, details of the computational method are presented. 
The distribution of free carriers and the band structures 
of graphene layers under different gate voltages are pre¬ 
sented in Sec. III. The effective dielectric constant of 
the /i-BN thin barrier in FETT is analyzed in Sec. IV. 
Finally, we give a summary on our results in Sec. V. 


II. SIMULATION METHOD 

The structure of the FETT used in experiments® is 
schematically shown in Fig. 1(a). Doped silicon is used as 
the gate electrode, which is separated from the current- 
carrying GrxI/i-BNjGrB trilayer channel by insulating sil¬ 
icon oxide (~ 300 nin) and /i-BN thin films (~ 50 nm). 
The gate voltage is applied between the doped silicon and 
the trilayer. 

In our DFT simulations, the doped silicon gate elec¬ 
trode was replaced by a semi-infinite ideally metallic 
medium for 2 < —zi [Fig. 1(b); the 2 -direction is perpen¬ 
dicular to the GrT|/i-BN|GrB trilayer], and the dielectrics 
between the gate electrode and the GrxI/i-BNjGrB tri¬ 
layer were replaced by vacuum to save computational re¬ 
sources. Thus the Grx|/i-BN|GrB trilayer is sandwiched 
by a semi-infinite vacuum medium (2 > 21 ) and an ide¬ 
ally metallic medium (2 < — 21 ), see Fig. 1(b). Note that 
the semi-infinite vacuum and metallic media are not ex¬ 
plicitly included in calculations, instead they are used as 
boundary conditions for the Hartree potential (I4r), 

dVn 

= 0 , = 0 . ( 1 ) 

We adopted the ESM inethod^^ to solve the Hartree po¬ 
tential. The purpose to employ this method is twofold: 
(i) Long-range Goulomb interactions with periodic im¬ 
ages are avoided, (ii) The non-periodic boundary condi¬ 
tion of Eq. (1) for the Hartree potential can be imposed. 

A gate voltage can be simulated by adding extra elec¬ 
trons or holes to the Grx|/i-BN|GrB trilayer;^^ the areal 
density of free carriers is proportional to the gate volt¬ 
age. In experiments the carrier density in graphene can 
be tuned by gate voltage up to 10^®cm“^, which is 
equivalent to 5.24 x 10“® electrons per primitive unit 
cell of graphene, or an electric displacement field of 
Dg = 0.016 C/m^ in the dielectrics. 

The measured interlayer distance between /i-BN and 
graphene layers^® is 3.32 A, and the interlayer distance 
between /i-BN layers 3.33 A is taken as half of the lattice 


constant c of crystalline /i-BN.^^ In-plane lattice con¬ 
stants for both graphene and /i-BN were set to 2.46 A 
due to the small lattice mismatch between them. In 
the x-y plane a periodic boundary condition was ap¬ 
plied with a dense 155 x 155 Monkhorst-Pack^® fc-mesh. 
The cutoff energy for wave functions and the Methfessel- 
Paxton^® spreading energy were taken to be 70 Ryd and 
10“®Ryd, respectively. The Perdew-Burke-Ernzerhof 
parameterization^^ of the generalized gradient approx¬ 
imation (GGA) to the exchange-correlation functional 
and Trouiller-Martins norm-conserving pseudopotentials 
were used. DFT calculations were performed using the 
Quantum ESPRESSO package^®. 

Charge density in the Gi'Bj/i-BNjGrx trilayer was cal¬ 
culated self-consistently for different gate voltages. In 
order to illustrate the distribution of free carriers across 
the trilayer, it is convenient to integrate the charge den¬ 
sity in the x-y plane and define a charge density along 
2 -direction p{z). The density of free carriers is defined as 
charge density difference for a device under a finite gate 
voltage Vg with respect to that under zero gate voltage. 
Compared to gate voltages, the electric displacement field 
Dg in the dielectrics is a more convenient quantity, be¬ 
cause it is independent of the thickness or the dielectric 
constant of the dielectrics. Gate voltages are expressed 
as Dg hereafter. 
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FIG. 2. (Color Online) (al) The distribution of free carriers 
under Dg = ±3 x 10~® C/m^ and (a2) the charge density 
under zero gate voltage of a Grx [monolayer /i-BN|GrB system, 
(b) The free carrier distribution in a graphenejmultilayer h- 
BNj graphene FETT under Dg = 3 x 10~® C/w? with the 
positions of the top and (c) the bottom graphene layer aligned, 
respectively. 
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III. RESULTS 


A FETT with a monolayer /i-BN barrier is used as an 
example to illustrate the charge and free carrier densi¬ 
ties in such devices.. The charge and free carrier den¬ 
sity on each layer can be analyzed using the Bader 
decomposition/® in which boundaries between atoms are 
defined by zero flux surfaces. There are three peaks 
in the charge density curve p{z), corresponding to the 
Grx, /i-BN, and Gre layers, respectively; cf. Fig. 2(a2). 
Free carrier densities Ap{z), defined as the change in 
charge density induced by gate voltages, under Dg = 
±0.003 G/m^ are shown in Fig. 2(al), where it can be 
seen that they have opposite sign but almost the same 
amplitude, in accordance with the electron-hole symme¬ 
try of graphene in the vicinity of the Fermi energy. For 
the one-dimensional charge density p{z), atomic layers 
are divided by minima of p[z), as denoted by vertical 
dashed lines in Fig. 2(a2) (zero flux surfaces shrink to 
points for the one-dimensional charge density.) A large 
portion of the free carriers accumulate on Gre at the side 
closer to the gate electrode; thus the free carrier density 
on GrB is larger than that on Grx. The asymmetric 
shape of Ap{z) with respect to the central /i-BN layer 
is a consequence of the asymmetric environment of the 
GrB|/i-BN|GrT trilayer, with the metallic gate electrode 
on one side and vacuum on the other side. 

The h-BN barrier can screen the applied gate voltage 
by developing an internal electric polarization, which re¬ 
duces the effect of the gate voltage and the free carrier 
density in Grx. We performed calculations on devices 
with a h-BN barrier thickness of up to eight monolayers. 
The free carrier density for these devices as a function 
of z for Dg = 3 X 10“^ G/m® is shown in Figs. 2(b) and 
2(c), shifted to make the positions of the Gi't (b) or 
GrB (c) layers coincide. The distribution of free carri¬ 
ers on Grs is almost independent of the thickness of the 



FIG. 3. (a) Band structure of a FETT with five layers of 
h-BN as barrier at Dg — 0.016 C/m^, and (b) the difference 
in chemical potential between Gre and Grx as a function of 
gate voltage (the dashed line is a guide to the eye). 


h-BN barrier, cf. Fig. 2(c). The amplitude of the elec¬ 
tric polarization in the h-BN barrier shrinks for thicker 
barriers (Fig. 2(c)), and so does the free carrier density 
in Grx (Fig. 2(b)), indicating that thicker h-BN barriers 
provide stronger screening. 

The chemical potential of graphene, defined as the 
energy difference between the Fermi energy and the 
charge neutrality point, can be efficiently tuned by using 
gate voltages because of the small DOS near the Fermi 
energy.®^® The band structure of a GrxI/i-BNjGrB trilayer 
in the vicinity of the Fermi energy has contributions only 
from graphene layers, because /i-BN is a wide gap insu¬ 
lator with a calculated energy band gap of 4.77 eV. The 
band structure of graphene shows a tiny energy gap of 
0.05 eV [Fig. 3(a)] induced by interaction with /i-BN.®® A 
FETT with a h-BN barrier thicker than one monolayer 
shows no hybridization between Gcb and Grx in its band 
structure. As an example, the band structure of a FETT 
with a five-layer-thick h-BN barrier for Dg = 0.016 C/m® 
is shown in Fig. 3(a), where bands from Grx and GrB 
shift upward rigidly by ^ 0.1 eV and ~ 0.2 eV with re¬ 
spect to the Fermi energy, respectively. The bands from 
Grx are always shifted away from the Fermi level by a 
larger energy than those from GrB [Fig. 3(b)]. 


IV. DISCUSSION 


Macroscopic electrostatic models have been used to 
calculate the electrostatic potential and free carrier den¬ 
sity in graphene based FETT.®’®^^®® In these models, the 
h-BN barrier between graphene layers is assumed to have 
a dielectric constant equal to that of the h-BN crystal. 
However, interfaces can exhibit significantly different di¬ 
electric properties compared to the bulk,^^ and the di¬ 
electric properties of a few-layer-thick h-BN barrier sand¬ 
wiched between two graphene layers need to be revisited. 
We seek to find the (effective) dielectric constant cbn of 
a h-BN barrier in a FETT and compare to its bulk value, 
in order to investigate if the dielectric constant is still a 
valid physical concept for thin h-BN films and to inves¬ 
tigate whether cbn of a h-BN barrier is modified by the 
interfaces with graphene layers. 

The dielectric constants of h-BN thin films (ebn) sand¬ 
wiched between graphene layers can be deduced from the 
electric field inside the /i-BN thin film (Ejnside) and the 
electric polarization (P) induced by gate voltages. 


ebn 


Ainside ± .H/cq 
.^inside 


( 2 ) 


where Emside is equal to the slope of the self-consistent 
Kohn-Sham effective potential in the z-direction. Einside 
is also related to the difference between chemical poten¬ 
tials of graphene layers Ap by ePinsided = Ap, where d 
is the distance between the two graphene layers. 

The polarization can be expressed as the summation 
of centers of Wannier functions according to the mod¬ 
ern theory of polarization.®^’®® In practice we followed 
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the method proposed in Refs. 26 and 27 which extends 
the Wannier function theory of polarization^'*’^®’^®’^® to 
metal-insulator heterostructures. 

The hybrid Wannier functions^®’^^, which are exponen¬ 
tially localized in the z direction but Bloch-like in the 
x-y plane, were constructed using the parallel-transport 
method.The first Brillouin zone was sampled by dis¬ 
crete fc-points of the type k = +j b, where the vectors 
kj_ form a 7V_l x N± uniform mesh in the x-y plane, and 
b = (0, 0, 27r7V||/L) is along the z-direction with L the 
height of the unit cell and iV|| the number of fc-points 
along the z-direction. The total number of fc-points is 
We used lVj_ = 7 and iV|| = 3, which are suffi¬ 
cient to converge the polarization. 

The matrices Mmn(k) = (Mm,k|u„,k+b) were con¬ 
structed where |um.k) is periodic under lattice trans¬ 
lations and n is the band index; |um,k) is normalized 
such that (um k) = 1. Singular value decomposi¬ 
tion of each M matrix was done utilizing the LAPACK 
library: M = C/EPl where U and V are complex uni¬ 
tary matrices and E is a diagonal matrix with diago¬ 
nal elements very close to 1 for small b. A new matrix 
M = UV^ was constructed corresponding to each M ma¬ 
trix, and a global matrix A(kj.) = rijio M(k± jh) 
was then constructed for each kj^ point. The centers 
of hybrid Wannier functions were calculated as Zm = 
(—L/27r) Im[ln Am], where the Am are the eigenvalues of 
A. 

The number of bands considered to construct the hy¬ 
brid Wannier functions, i.e., the dimension of the M ma¬ 
trices, is equal to four for each graphene or /i-BN atomic 
layer, the number of occupied bands in most of the first 


Brillouin zone except for the small portion near the K- 
point (see Fig. 3). This choice does not affect the cal¬ 
culations of the polarization inside /i-BN thin films be¬ 
cause the bands near the Fermi energy are contributed 
by graphene layers. 

Four of the resulting hybrid Wannier functions can be 
assigned to each graphene or /i-BN layer according to the 
positions of their centers. Two of them are very close to 
the atomic plane (within 10“^ A) and the other two are 
located about 0.4 A above and beneath the atomic plane 
respectively. The center of charge for each atomic layer 
is equal to the average value of the corresponding four 
Wannier functions. The dipole moment corresponding 
to each /i-BN layer was calculated using the shift of the 
center of charge under an electric field, and the polariza¬ 
tion was calculated with the thickness of each /i-BN layer 
set to be 3.33 A. 

The calculated polarization of the /i-BN layer adjacent 
to the graphene layers is almost the same as /i-BN layers 
deeply inside [see Fig. 4(a)], indicating that interface with 
graphene layers has little effect on the dielectric proper¬ 
ties of /i-BN thin films. As a result, the average dielectric 
constant for /i-BN thin films embedded by graphene lay¬ 
ers is independent of the thickness; see Fig. 4(b). 

The effect of an interface with graphene on the dielec¬ 
tric properties of /i-BN can be analyzed using the inter¬ 
face charge redistribution, denoted as Apintf and defined 
as the difference in charge redistribution at the interface 
with respect to that in the bulk (denoted as Apbuik)- In 
practice we used the charge redistribution at the center 
of the /i-BN as the bulk, as shown in Fig. 5. Thus one 
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FIG. 4. (a) The polarization of /i-BN layers in a FETT with 
a five-layer /i-BN thin film induced by an gate voltage of 
Dg — 6.114 X 10“®C/m^. (b) The calculated average rela¬ 
tive dielectric permittivity of /i-BN thin films embedded by 
graphene layers as a function of the number of /i-BN atomic 
layers, and the dashed line denotes the relative dielectric per¬ 
mittivity of bulk /i-BN. 
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FIG. 5. (Color online) Interface charge redistribution Apintf 
in a FETT device with eight-layer /i-BN for Dg = 0.016 C/m^, 
where Apintf is the total charge redistribution Ap minus 
the bulk charge redistribution Apbuik at the center of /i-BN; 
Apbuik is denoted by the patterned rectangle in the top panel. 
The absolute value of Apintf is plotted on a logarithmic scale 
in the bottom panel. The positions of the atomic planes of 
Gru, Grx and the two interface /i-BN layers BNi and BNs 
are denoted by vertical dashed lines. 
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can obtain Apintf by subtracting Apbuik from the total 
charge redistribution. Because /i-BN is a wide-band-gap 
insulator, Apintf should decay exponentially away from 
the interface. The strength of the interface effect is de¬ 
termined by the amplitude of Apintf near the interface. 

An example the charge redistribution in a FETT de¬ 
vice with eight-layer /i-BN for Dg = 0.016 C/m^ is shown 
in Fig. 5. The difference Apintf is large near the graphene 
layers, decays quickly into /i-BN, and becomes invisible 
after crossing the first /i-BN atomic layers (BNi and BNg 
in Fig. 5). The amplitude of Apintf is also presented on a 
logarithmic scale in the lower panel of Fig. 5. The decay 
of Apintf into /i-BN is approximately exponential. Most 
importantly, the amplitude of Apintf is about 50 times 
smaller than that of Apbuik between BNi and BNg in 
Fig. 5: Apintf is very small inside the atomic plane of the 
first /i-BN layer. As a result, any interface effect on the 
dielectric properties of /i-BN is very weak, which explains 
why the calculated dielectric constant of /i-BN in FETT 
devices is close to the bulk value. 
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FIG. 6. (Color online) (Upper panel) The atomic strncture 
of a hydrogen satnrated Si(lll) thin film with a thickness 
of abont 2 nm. Si and H atoms are represented by large 
blue and small pink spheres, respectively. Solid lines denote 
the boundary of the nnit cell. (Lower panel) The shift of 
the hybrid Wannier functions indnced by an electric field of 
0.0385 V/A along the z-direction. Red circles (o) and blue 
discs (•) represent hybrid Wannier fnnctions with higher and 
lower polarizability, respectively. 


We plotted in Fig. 6 the shift of each hybrid Wan¬ 
nier functions induced by an electric field of 0.0385 V/A 
along the z-direction, where the electric field was applied 
using the ESM method using a metal | slab | metal config¬ 
uration. The hybrid Wannier functions can be divided 
into two categories according to their polarizability. The 
hybrid Wannier functions located at canted Si-Si bonds 
with respect to the z-direction (denoted as • in Fig. 6) 
exhibit lower polarizability, and they show negligible de¬ 
viations at the surface. On the other hand, the hybrid 
Wannier functions located at parallel Si-Si or Si-H bonds 
with respect to the z-direction (denoted as o in Fig. 6) 
exhibit higher polarizability. We also observed that the 
hybrid Wannier functions of Si-H bonds at the surface 
show a polarizability 12% lower than those correspond¬ 
ing to parallel Si-Si bonds. The 12% lower polarizability 
of the Si(lll) surfaces shown in Fig. 6 is not as severe 
as reported by previous studies.In those studies, the 
macroscopic polarization and dielectric constant were ob¬ 
tained after a smoothing procedure. The purpose of the 
smoothing procedure is to eliminate the dielectric nonlo¬ 
cality, but this procedure reduces the dielectric permit¬ 
tivity of surfaces artificially because the dielectric con¬ 
stant of vacuum is smeared into the surface. 


V. SUMMARY 


The distribution of free carriers and the band struc¬ 
ture of graphene layers in graphene based FETT have 
been simulated using the DFT-I-ESM method. The di¬ 
electric properties of /i-BN thin films sandwiched between 
graphene layers in EETT were investigated using the the¬ 
ory of microscopic permittivity and found to have a di¬ 
electric permittivity close to that of crystalline /i-BN. The 
small amplitude of interface charge redistribution inside 
the atomic plane of the first /i-BN layer proves that the ef¬ 
fect of the interface with graphene on the dielectric prop¬ 
erties of /i-BN is weak. 

In this study we have demonstrated the DFT-I-ESM 
method as a promising approach to simulate field-effect 
devices with a planar structure. Once the charge density 
and effective potential of a field-effect device are self- 
consistently obtained, the scattering of transport elec¬ 
trons and electric conductivity can be calculated using 
scattering theory. 


We also compared the /i-BN thin films with silicon thin 
films, because the latter are known to exhibit lower di¬ 
electric permittivity than for the corresponding bulk.^^’^® 
The in-plane lattice constant of Si(lll) slabs was chosen 
to be the experimental lattice constant of /cc-Si. Dan¬ 
gling bonds on both surfaces are saturated by hydrogen, 
as shown in the upper panel of Fig. 6, and the Si-H bond 
length is 1.50 A after structure optimization. 
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